Self-organized critical earthquake model with moving boundary 
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A globally driven self-organized critical model of earthquakes with conservative dynamics has 
been studied. An open but moving boundary condition has been used so that the origin (epicenter) 
of every avalanche (earthquake) is at the center of the boundary. As a result, all avalanches grow in 
equivalent conditions and the avalanche size distribution obeys finite size scaling excellent. Though 
the recurrence time distribution of the time series of avalanche sizes obeys well both the scaling 
forms recently observed in analysis of the real data of earthquakes, it is found that the scaling 
function decays only exponentially in contrast to a generalized gamma distribution observed in the 
real data analysis. The non-conservative version of the model shows periodicity even with open 
boundary. 

PACS numbers: 05.65.+b 91.30.Dk 64.60.Ht 89.75.Da 



Because of the devastating effects of earthquakes on 
human life and wealth, understanding the properties, 
behavior and statistics of earthquakes as well as their 
predictions continue to remain a challenge to scientists. 
Over a long time attempts have been made to explain the 
earthquake dynamics as a scale invariant process. For ex- 
ample, Gutenberg-Richter distribution law for the earth- 
quake magnitudes Omori's law for the frequencies 
of after shocks 0] as well as recent analysis of recur- 
rence time distributions [H, 0, 0, El , fractal distribution 
of epicenters 0, S] , power law distribution of the spatial 
distances between epicenters of successive earthquakes 
, and associating a scale-free network with the tempo- 
ral behaviour of earthquakes [lo| . all support the view 
point that earthquakes are indeed scale invariant. On 
the other hand theoretically, the well known Burridge- 
Knopoff (BK) model views the slow creeping of the con- 
tinental plates along the fault lines as a stick-slip pro- 
cess [llj . About two decades ago, Bak et. al. while 
introducing the idea of Self-Organized Criticality (SOC) 
suggested that the phenomenon of earthquakes may be 
looked upon as a SOC process since there is nobody to 
control the nature to generate long r ang e spatio-temporal 



correlations or scale-invariance 

In this paper we study a SOC model of earthquakes 
and present numerical evidence to argue that within the 
frame-work of this model the earthquake dynamics is in- 
deed scale-invariant. In particular, we show that the two 
recently used scaling procedures for analyzing the real 
data of earthquakes work well for our model. 

The well known Gutenberg-Richter (GR) law says that 
the number of earthquakes N(m) of magnitude at least 
m decays exponentially with m as: 



log N(m) = C\ — C2m. 



(1) 



On the other hand magnitude of an earthquake varies 
logarithmically with the amount of energy released: 
log E(m) — C3 + c/jn. Eliminating m one gets, logiV = 
c i — (c2/c4)logE + (c2C3)/c4. This implies that the cu- 
mulative number N(E) of earthquakes of energy at least 



E decays like a power law as: 

N(E) oc E~ b 



(2) 



where b = C2/C4. Therefore the probability density of 
earthquakes varies as: Prob(.E') ex dN(E)/dE oc E~ x ~ h . 
Another important empirical observation is the Omori 
law which states that the frequency of after shocks decays 
with time as a power law: D(t) oc t~ 7 . 

Let the earthquakes be measured with an accuracy m c 
so that all earthquakes of magnitude greater than m c = 
log 10 (s c ) are detected and let them be ordered in a time 
sequence so that the i-th earthquake occur at the time 
ti. The recurrence time is then denned as Ti = ti — tj_j.. 
Bak et. al. analyzed the real data of the earthquakes 
occured in Southern California by dividing this region 
into a grid of cell size L degrees. Considering main events, 
after shocks and fore shocks on the same footing Bak, 
Christensen, Danon and Scanlon (BCDS) claimed that 
the recurrence time follows an universal scaling function 



3] 



Prob(r, L, s c ) ~ t 



(3) 



where 6, 7 and df are the GR exponent, Omori exponent 
and the fractal dimension of the distribution of epicenters 
and J-{x) is an universal scaling function. The scaling 
factor s\/L ds is the mean recurrence time for the earth- 
quakes having sizes at least s c which originated from a 
cell of size L. On the other hand Corral used a single 
parameter R for scaling, which is the rate of occurrence 
of the earthquakes [5|, l6| : 



Prob(T, R) - RG{Rt) 



(4) 



where G(x) is another universal scaling function having 
the form of a generalized Gamma function. 

Bak and Tang devised a SOC model of earthquakes 
by studying a simpler version of the two-dimensional BK 
model |13| . The essential simplification is in treating the 
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FIG. 1: Four examples show the positions of the avalanche 
origins (shaded circle) and the corresponding boundaries sites 
(filled circles) on a 8 x 8 square lattice. 




accumulated local force as a scalar as well as consider- 
ing the two-dimensional system of blocks located at fixed 
positions at the sites of a regular lattice like a discrete 
space-time but continuous spin cellular automaton. 

Olami, Feder and Christensen (OFC) studied the non- 
conservative version of the SOC model of earthquakes 
14j . Every site of a square lattice is assigned a continuous 
variable / representing the accumulated local force at 
that site. The system is globally driven, implying that 
in the inactive state of no avalanches (earthquakes) the 
forces at all sites increase steadily with an uniform rate. 
A threshold value f c of the forces exists for the stability 
of all sites. A site relaxes with probability one when 
fi.j > fc- In a relaxation the force at the site is reset 
to zero and a fraction of the force is transmitted to each 
neighbor: 



If, fi,j > fc, then fij -> and 

fi±l,j±l — > fi±l,j±l + OtfiJ. 



(5) 



Consequently, forces at some of the neighbors may exceed 
the threshold and in turn they also relax - thus a cascade 
of site relaxations propagates in the system, causing an 
avalanche. The parameter a varies continuously within 
the range < a < 1/4 Q. The dynamics in OFC 
model is conservative for a = 1/4 and non-conservative 
of a < 1/4. A critical value of a c such that the system 
is in a sub-critical state for a < a c and in a critical state 
for a > a c has been suggested for a c « 0.05 around 
0.20 0, = l/4 0anda multifractal scaling in [l7|- 

We argue that assigning a fixed boundary in the SOC 
models of earthquakes is rather artificial. In nature there 
is no fixed boundary for the earthquakes which absorbs 
earth's vibrations, the seismic waves propagate in all 
directions till they slowly damp out at long distances. 
Presence of a fixed boundary introduces a strong non- 



FIG. 2: (Color online) The avalanche size distributions for 
three different system sizes L — 64 (black), 128 (red) and 256 
(blue) have been plotted on a double logarithmic scale in (a) . 
The finite size scaling of the same data is shown in (b). 



uniformity in the system i.e., all measurable quantities 
show strong dependence on the distance from the bound- 
ary. This effect is present in both conservative as well as 
non-conservative versions of the OFC model, but it is so 
strong in the latter case that even arriving at the sta- 
tionary state becomes very difficult T5J. It is therefore 
desirable that all avalanches are on the same footing with 
respect to the boundary and at the same time the origin 
of the avalanche should be at the farthest interior point 
of the system. 

This argument prompted us to formulate a new bound- 
ary condition. Here, a globally fixed set of lattice sites 
does not constitute the boundary for all avalanches. In 
contrast, boundaries are different for different avalanches 
depending on the positions of the avalanche origins, and 
its position is constantly moved from one avalanche to 
the other. 

First we make the square lattice periodic in both di- 
rections to get the topology of a torus. An arbitrary 
random distribution of forces fij, drawn from a set of 
independent and identically distributed random numbers 
within {0, 1} are assigned at all L 2 sites. The maximum 
force fmax among all L 2 sites is found to be at some spe- 
cific location (i , j a ) and the difference from the threshold 
force is estimated: 6 = f c — fmax- Forces at all sites are 
then enhanced by 8 so that at the origin (i ,jo) the force 
reaches the threshold f c . The avalanche then starts from 
the origin and a cascade of relaxations propagates away 
from the origin. 
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Now, for this avalanche, we select a specific set of 
lattice sites as the boundary such that the origin is at 
the center position with respect to these boundary sites. 
More precisely, on a L x L square lattice and with respect 
to the origin located at (i ,j ) the boundary sites form 
two transverse rings on the torus defined by one column 
and one row of lattice sites as (Fig. 1): 

i = i + L/2 mod(L) and 

3 = Jo + L/2 mod(L). (6) 

When a site adjacent to the boundary relaxes, it transfers 
afij force to every non-boundary neighbor but no force 
to the neigbor on the boundary. Therefore corrsepond- 
ing to each boundary neighbor afij disappears from the 
system and in this way the system looses force. 

Since the system is otherwise periodic in all direc- 
tions all lattice sites are equivalent. Consequently all 
avalanches are also equivalent since all of them grow in 
similar surroundings. In a way this is similar to elim- 
ination of surface effects in a finite size system. Sur- 
face profiles for the averaged force per site (/), number 
of avalanche origins at each site (e) and average size of 
the avalanche per site (s) show uniform flat surfaces but 
within a very small fluctuation for all sites within the lat- 
tice L x L. We are also studying other numerically chal- 
lenging problems of SOC using moving boundary condi- 
tion. 

Since in a single relaxation, the force at the site is re- 
duced to zero, it creates the possibility that more than 
one site (typically two) can reach the threshold simulta- 
neously. However, such situations occur with very low 
probability and in these cases we choose randomly one 
of the sites as the origin and construct boundaries with 
respect to this site but relaxation starts from both the 
unstable sites. Since the forces are continuously varying 
real numbers, t he p recision of the numbers is important 
as observed in [18|. To ensure that the system has in- 
deed reached the stationary state we calculated the aver- 
age avalanche size (s(L)) for every 10000 avalanches and 
monitored its variation with time. This quantity first 
grows with time but eventually saturates. Repeating this 
calculation for different system sizes it is observed that 
the relaxation time grows as L 2 . 

First we present the results for the conservative case 
of a = 1/4. The avalanche size s is the total number of 
relaxations in an avalanche and represents the total en- 
ergy release in our model earthquake. Prob(s, L) is the 
probability that a randomly selected avalanche has size 
s. While for the infinitely large system size the distribu- 
tion should indeed be a simple power law, for the finite 
size systems, a finite size scaling of the distribution is 
required: 

Prob(s,L) ~ L-»H(s/L v ) (7) 

where the scaling function Tt(x) ~ x^ 1 ^ for x — > and 
for x » 1, TL(x) decreases faster than a power law so 




t(lV) 



FIG. 3: (Color online) Scaling of the waiting time distribution 
by the (a) Corral method (b) BCDS method. Symbols used 
for: L = 64 (circle), 128 (square), 256 (triangle) and for s c = 
(black), 8 (red), 64 (blue) and 512 (magenta). Values of the 
scaling exponents used in (b) are df =1.67 and 6=0.29. The 
continuous line is the best fit by the functional form in Eqn. 
(8). 
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FIG. 4: (Color online) The average recurrence time (r(L, s c )} 
has been plotted for different values of s c and multiplied by 
the system size dependent factor L d f . On increasing system 
size the plot approaches to the variation mentioned in Eqn. 
(10) with b = 0.30. 

that, b — fi/v — 1. The system size dependence of the 
average avalanche size and durations are observed to be 
(s(L)) ~ L 226 and (T(L)) ~ L 63 . This shows that 
the avalanche dynamics is sub-diffusive. We believe that 
this is due to fact that force is always reset to zero in 
a relaxation which initiates more relaxations and thus 
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increases the size of the avalanche. 

In Fig. 2(a) we show the plot of avalanche size distri- 
bution for three different system sizes L — 64, 128 and 
256 on the double logarithmic scale. All of them have 
very large portions of straight regions starting from very 
small sizes to the cut-off sizes. A scaling of the data with 
an excellent data collapse is shown in Fig. 2(b) yielding 
the values of v = 3.02 and fi — 3.78 giving b 0.26. Such 
a good power law behaviour as well as the excellent finite 
size scaling have been achieved only due to the moving 
boundary condition where all lattice sites as well as the 
avalanches are equivalent and have not been observed in 
fixed boundary cellular automata models of earthquakes 
before 0, EH. 

Since we have assumed that forces at all sites in- 
crease uniformly at unit rate, the time difference be- 
tween successive avalanches is exactly equal to S. With 
this definition, the recurrence time distribution (RTD) 
Prob(r, s c , L) has been calculated for different system 
sizes L as well as different s c values. The effects of s c and 
L on RTD are competitive. For s c = 0, the RTD is sim- 
ply the distribution of force increments S only. Since the 
probability of occurrence of an avalanche of size at least 
s c decreases with s c , for a fixed L the recurrence time in- 
creases with increasing s c . On the other hand for a fixed 
s c , since the maximum of the avalanche sizes increases 
with L, the probability of occurrence of an avalanche of 
size at least s c increases with increasing L. Consequently 
the recurrence time decreases with increasing L. 

In Fig. 3(a) we show an unified scaling of twelve differ- 
ent plots with the minimal value of the avalanche sizes 
measured s c = 0, 8, 64 and 512 for three different sys- 
tem sizes L =64, 128 and 256. Logarithmic binning is 
used for coarse-graining of the data. The average wait- 
ing time (t(L, s c )) is calculated for each plot. Following 
Eqn. (4) we then scale every plot with corresponding 
(t(L,s c )) — 1/R and observe an excellent collapse of 
all twelve plots. This confirms the Corral scaling in our 
model. We tried to verify the Corral scaling form: 

Q(x) ~ x ~ ai exp{-a 2 x a3 ) (8) 

and obtained ai = 0.003, a 2 = 1-02 and a 3 = 0.99 com- 
pared to a\ — 0.33, a 2 = 0.63 and = 0.98 observed in 
[5|. The exponential tail in 0(x) is consistent with the 
Gamma distribution observed by Corral but the observed 
power law decay component for small values of waiting 
times is rather absent in our model. 

To see if BCDS scaling is also valid for our model, we 
plotted Prob(r, L, s c ){s b c /L d f) vs. rL d f/s b c and obtained 
a scaling form like: 

Prob(r,L, Sc )^-~ ?i(t^). (9) 

Here also we see a very good collapse of the nine sets 
of data for three system sizes L = 64, 128 and 256 and 



for s c = 8, 64 and 512. The scaling exponents that gave 
the best collapse were tuned to df — 1.67 and b = 0.29. 
The best fit with the functional form in Eqn. (8) gives 
ai = 0.001, a 2 = 3.21 and = 0.99 again showing an 
exponential tail similar to that obtained from real data 
analysis [B[ but without any power law component. 

We therefore conclude that both the scaling forms used 
by Corral as well as BCDS are valid for scaling of the 
RTD data in our model. The scaling functions in both 
cases were observed to be very close to simple exponential 
decay and the power law part representing the RTD for 
small values of the recurrence times turned out to be 
absent. This result may also be compared with two recent 
analytical calculations: (i) a pure exponential decay of 
the RTD [3] (ii) an approximate unified law compatible 
with the empirical observations incorporating the Omori 
law [U. 

For Corral's analysis it is the single parameter scaling 
i.e., the mean recurrence time (r(L, s c )}. However this 
parameter in turn also depends jointly on the another 
two competitive parameters of the distribution, namely 
the system size L and the avalanche size cut-off s c in the 
following way: 

(t(L,s c ))<x^j. (10) 

To check if it is really true we plotted (t(L, s c ))L df with 
respect to s b c for L — 32, 64, 128 and 256 using df = 1.67 
in Fig. 4. A nice collapse of the data for the four differ- 
ent system sizes are observed for small and intermediate 
values of s c . Collapse of the data between two successive 
system sizes increased with the system size. The slope 
of the curve in the longest straight region corresponds to 
b= 0.30. 

Finally, we studied the OFC model using values of 
a < 1/4 again on a square lattice of size L using open 
but moving boundary condition. To our surprise we see 
that the dynamics become periodic after a short relax- 
ation time of the order of L 2 . This is checked by looking 
at the 'hamming distance'. Starting from a random dis- 
tribution of forces as before we skip some 10L 2 initial 
avalanches and store the force configuration in an array 
f store- After that at the end of every avalanche we cal- 
culated the maximal site difference max\fij — f s tore(hj)\ 
and measure the time when this maximal difference be- 
comes less than a small number e = 10~ 12 . The periodic 
time is of the order of L 2 but less than it, and found to 
depend on the initial distribution of force values. 

To summarize, we have studied in a model the scale 
invariance properties observed in the real data of earth- 
quakes over last several years by different groups. More 
specifically we studied a self-organized critical model of 
earthquakes using a square lattice cellular automaton. 
Using a moving boundary condition we could eliminate 
all boundary effects. We first observe that the avalanche 
size distribution of this model follow excellent finite size 
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scaling. Further, the recurrence time distribution was 
analyzed in two ways, i.e., using Corral as well as BCDS 
scalings. We observe that our simulated data of the RTD 
support both scalings very well which leads us to con- 
clude that the mean recurrence time is actually a joint 
function of both the system size as well as the avalanche 
size cut-off as used to measure the waiting times. 
Electronic Address: manna@boson.bose. res. in 
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